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Abstract. Recent analyses of the helioseismic observa- 
tions indicate that the previously observed surface tor- 
sional oscillations with periods of about 11 years extend 
significantly downwards into the solar convective zone. 
Furthermore, there are indications that the dynamical 
regimes at the base of the convection zone are different 
from those observed at the top, having either significantly 
shorter periods or non-periodic behaviour. 

We propose that this behaviour can be explained by 
the occurrence of spatiotemporal fragmentation, a crucial 
feature of which is that such behaviour can be explained 
solely through nonlinear spatiotemporal dynamics, with- 
out requiring separate mechanisms with different time 
scales at different depths. 

We find evidence for this mechanism in the context 
of a two dimensional axisymmetric mean field dynamo 
model operating in a spherical shell, with a semi-open 
outer boundary condition, in which the only nonlinearity 
is the action of the azimuthal component of the Lorentz 
force of the dynamo generated magnetic field on the solar 
angular velocity. 
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1. Introduction 

Recent analyses of the helioseismic data, from both the 
Michelson Doppler Imager (MDI) instrument on board the 
SOHO spacecraft (Howe et al. 2000a) and the Global Os- 
cillation Network Group (GONG) project (Antia & Basu 
2000) have provided strong evidence which indicates that 
the earlier observed time variations of the differential ro- 
tation on the solar surface, the so called 'torsional oscil- 
lations' with periods of about 11 years (e.g. Howard & 
LaBontc 1980; Snodgrass, Howard & Webster 1985; Koso- 
vichev & Schou 1997; Schou et al. 1998), penetrate into 
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the convection zone, to depths of at least 8 percent in 
radius. 

Furthermore, these data have provided some evidence 
to suggest that variations in the differential rotation are 
also present around the tachocline at the bottom of the 
convection zone (Howe et al. 2000b) . An important feature 
that distinguishes these variations from those observed at 
the upper parts of the convection zone is that they possess 
markedly different modes of behaviour: either possessing 
distinctly lower periods (of ~1.3 years), or being non- 
periodic (Antia & Basu 2000). Clearly, to firmly establish 
the precise nature of these variations, future observations 
are required. Whatever the outcome of such observations, 
however, both these sets of results point to the very in- 
teresting possibility that the variations in the differential 
rotation can have different periodicities/behaviours at dif- 
ferent depths in the solar convection zone. It is therefore 
important to ask whether such different variations can in 
principle occur at different parts of the convection zone 
and, if so, what could be the possible mechanism(s) for 
their production. 

The aim of this letter is to suggest that a natural 
mechanism for the production of such different dynami- 
cal modes of behaviour in the convection zone is through 
what we call spatiotemporal fragmentation, i.e. the occur- 
rence of dynamical regimes at (given) values of the control 
parameters of the system, which possess different tem- 
poral behaviours at different spatial locations. This is to 
be contrasted with the usual temporal bifurcations, with 
identical temporal behaviour at each spatial point, which 
require changes in parameters to occur. 

We find evidence for this mechanism in the context of 
a two dimensional axisymmetric mean field dynamo model 
in a spherical shell, with a semi-open outer boundary con- 
dition, in which the only nonlinearity is the action of the 
azimuthal component of the Lorentz force of the dynamo 
generated magnetic field on the solar angular velocity. The 
underlying angular velocity is chosen to be consistent with 
the most recent helioseismic data. 

In addition to producing different dynamical variations 
in the differential rotation, including different periods, at 
the top and the bottom of the convection zone, this model 
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is also capable of producing butterfly diagrams which are 
in qualitative agreement with the observations as well as 
displaying torsional oscillations that penetrate into the 
convection zone, as recently observed by Howe et al. 2000a 
and Antia & Basu 2000 and studied by Covas et al. (2000). 

2. The model 

We shall assume that the gross features of the large scale 
solar magnetic field can be described by a mean field dy- 
namo model, with the standard equation 



dB 



V x (u x B + aB - 77 V x B) 



(1) 



Here u = vcj> — 5V77, the term proportional to V77 repre- 
sents the effects of turbulent diamagnetism, and the ve- 
locity field is taken to be of the form v = vq + v', where 
vq = QorsmO, ilo is a prescribed underlying rotation law 
and the component v' satisfies 
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where D is the operator 

-p-s) and fiQ is the induction constant. The source of the 
sole nonlinearity in the dynamo equation is the feedback of 
the azimuthal component of the Lorentz force (Eq. (0)), 
which modifies only slightly the underlying imposed ro- 
tation law, but thus limits the magnetic fields at finite 
amplitude. The assumption of axisymmetry allows the 
field B to be split simply into toroidal and poloidal parts, 
B = Br + Bp = B<j) + V x A<f>, and Eq. @) then yields 
two scalar equations for A and B. Nondimcnsionalizing in 
terms of the solar radius R and time R 2 /rjo, where 770 is 
the maximum value of 77, and putting £1 = a = ctoa, 
V = VoV> B = BqB and v' = Q*Rv', results in a sys- 
tem of equations for A, B and v' . The dynamo parameters 
are the two magnetic Reynolds numbers R a = aoR/rjo 
and R u = fl*R 2 /i]o, and the turbulent Prandtl number 
Pr = vq/t}q. ST is the solar surface equatorial angular ve- 
locity and 77 = 77/770- Thus Uq and 770 are the turbulent 
magnetic diffusivity and viscosity respectively, i? w is fixed 
when 770 is determined (see Sect. ||), but the value of R a is 
more uncertain. The density p is assumed to be uniform. 

When attempting to model astrophysical systems, 
boundary conditions are often rather ill-determined. We 
try to make physically motivated choices. For our inner 
boundary conditions we chose B = 0, ensuring angular 
momentum conservation, and an overshoot-type condi- 
tion on Bp (cf. Moss & Brooke 2000). At the outer bound- 
ary, we used an open boundary condition dB/dr = 
on B and vacuum boundary conditions for Bp. The mo- 
tivation for this is that the surface boundary condition 
is ill-defined, and there is some evidence that the more 
usual B — condition may be inadequate. This issue has 
recently been discussed at length by Kitchatinov et al. 



(2000), who derive 'non-vacuum' boundary conditions on 
both B and Bp. 

Equations (|l|) and (0) were solved using the code de- 
scribed in Moss & Brooke (2000) (see also Covas et al. 
2000) together with the above boundary conditions, over 
the range ro < r < 1, < 9 < 7r. We set vq — 0.64; with 
the solar convection zone proper being thought to occupy 
the region r > 0.7, the region ro < r < 0.7 can be thought 
of as an overshoot region/tachoclinc. In the following sim- 
ulations we used a mesh resolution of 61 x 101 points, 
uniformly distributed in radius and latitude respectively. 

In this investigation, we took Slo to be given in 0.64 < 
r < 1 by an interpolation on the MDI data obtained 
from 1996 to 1999 (Howe et al. 2000a). For a we took 
a = a r (r)f(8), where f(8) = sin 2 8 cos 8 (cf. Rudiger &: 
Brandenburg 1995) and a r = 1 for 0.7 < r < 0.8 with 
cubic interpolation to zero at r — r$ and r = 1, with the 
convention that a r > and R a < 0. Also, in order to take 
into account the likely decrease in the turbulent diffusion 
coefficient 77 in the overshoot region, we allowed a simple 
linear decrease from 77 = 1 at r = 0.8 to 77 = 0.5 in r < 0.7. 

3. Results 

We calibrated our model so that near marginal excita- 
tion the cycle period was about 22 years. This determined 
Ru = 44000, corresponding to 770 ~ 3.4 x 10 11 cm 2 sec^ 1 , 
given the known values of $1* and R. The first solutions 
to be excited in the linear theory are limit cycles with 
odd (dipolar) parity with respect to the equator, with 
marginal dynamo number R a w —2.23. The even par- 
ity (quadrupolar) solutions are also excited at a similar 
marginal dynamo number of R a ~ —2.24. It is plausible 
that the turbulent Prandtl number be of order unity, and 
we set P r = 1. For the parameter range that we inves- 
tigated, the even parity solutions are nonlinearly stable. 
Given that the Sun is observed to be close to an odd (dipo- 
lar) parity state, and that previous experience shows that 
small changes in the physical model can cause a change 
between odd and even parities in the stable nonlinear so- 
lution, we chose to impose dipolar parity on our solutions. 

With these parameter values, we found that this 
model, with the underlying zero order angular velocity 
chosen to be consistent with the recent (MDI) helioseis- 
mic data, is capable of producing butterfly diagrams which 
are in qualitative agreement with the observations. An ex- 
ample is shown in Fig. [I]. (The polar branch is a little too 
strong, but this feature can be weakened by adjusting the 
latitudinal dependence of a (see also Covas et al. 2000).) 
The model can also successfully produce torsional oscilla- 
tions (see Fig. ||) that penetrate into the convection zone, 
similar to those deduced from recent helioseismic data 
(Howe et al. 2000a) and studied in Covas et al. (2000). 
We note, however, that an additional interesting feature 
of the present model is that the torsional oscillations have 
larger and more realistic amplitudes near the surface, of 
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the order of 1 nHz, much larger than was found previously 
using the boundary condition B — at the surface. 




Fig. 1. Butterfly diagram of the toroidal component of 
the magnetic field B at fractional radius r = 0.95. Dark 
and light shades correspond to positive and negative val- 
ues of 2?0 respectively. Parameter values are R a = —3.0, 
P r = 1.0 and R„ = 44000. 




Fig. 2. Variation of rotation rate with latitude and time 
from which a temporal average has been subtracted to re- 
veal the migrating banded zonal flows, taken at fractional 
radius r = 0.99. Darker and lighter regions represent posi- 
tive and negative deviations from the time averaged back- 
ground rotation rate. Parameter values are as in Fig. [I]. 

We found that the model is also capable of producing 
spatiotemporal fragmentation, near the base of the con- 
vection zone, hence resulting in oscillations in the differ- 
ential rotation with, for example, half the basic period. To 
demonstrate this, we have plotted in Figs. ||-|] the radial 
contours of the angular velocity residuals Sfl as a func- 
tion of time for a cut at latitude 30°, for several values 
of R a . As can be seen, for smaller values of R a (Fig. |J), 
we find torsional oscillations with the same period at the 
top and the bottom of the convection zone. As R a is in- 



creased (Figs. [I] and |5|), a spatiotemporal fragmentation 
occurs near the base of the convection zone, resulting in 
oscillations in the differential rotation with half the pe- 
riod of the oscillations near the top. For still higher values 
of R a , the temporal variations in the differential rotation 
at the base of the convection zone start to become non- 
periodic, which might be of relevance if the failure of Antia 
& Basu (2000) to find shorter period oscillations near the 
bottom of the convection zone should turn out to be cor- 
rect. We have also checked that the butterfly diagrams do 
not fragment and keep the same period independently of 
the depth and R a value, continuing to resemble Fig. |l]. 
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Fig. 3. Radial contours of the angular velocity residuals 
Sil as a function of time for a cut at latitude 30°. Param- 
eter values are R a = -3.0, P r = 1.0, R u = 44000. Note 
how the torsional oscillations are very coherent from top 
to the base of the dynamo region showing that only one 
period is present. 
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Fig. 4. Radial contours of the angular velocity residuals 
Sfl as in Fig. || for R a = —7.0. Note the emergence of 
spatiotemporal fragmentation towards the bottom of the 
convective zone, resulting in different periodicities there. 
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Fig. 5. Radial contours of the angular velocity residuals 
Sn as in Fig. §for R a = -10.0. 

This fragmentation is made more transparent in Fig. 
^ which shows the temporal oscillations in the angular 
velocity residuals Sfl at a fixed point, as R a is increased, 
illustrating the presence of period halving. 
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Fig. 6. 'Period halving' at r = 0.68 and latitude 30°. 
The panels correspond, from left to right, to R a values 
—3.0 and —7.0 respectively, and display increasing rela- 
tive amplitudes of the secondary oscillations. Remaining 
parameters values are as in Fig. ||. 



4. Discussion 

We have proposed spatiotemporal fragmentation as a nat- 
ural mechanism for producing different types of variation 
in the differential rotation at the top and the bottom of the 
convection zone. To demonstrate the occurrence of such 
behaviour, we have studied a solar dynamo model, with 
a semi-open outer boundary condition, calibrated to have 
the correct cycle period, with a mean rotation law given by 



recent helioseismic observations. We note in passing that 
in a few simulations performed with the boundary condi- 
tion B — at the surface, we have not so far found this 
phenomenon, although we cannot yet make a definitive 
statement on this point. In addition to producing butter- 
fly diagrams in qualitative agreement with those that are 
observed, as well as torsional oscillations that penetrate 
into the convection zone, we have shown that this model 
can also produce spatiotemporal fragmentation, resulting 
in different oscillatory modes of behaviour near the top 
and the bottom of the convection zone. 

We emphasize that the main aim of this letter is to 
propose a mechanism that can be expected to operate in 
general nonlinear dynamo settings, and which is capable 
of producing multiple periods and/or non-periodic oscilla- 
tions in parts of the convective zone. The specific results 
given here, such as the single period halving, are based 
on a particular dynamo model which inevitably includes 
many simplifying assumptions, not least of which is that 
the density is uniform. (It is unclear how the inclusion of 
a radial dependence p(r) would affect our results — we 
note that current solar dynamo models commonly take 
a uniform density.) We expect that the mechanism is of 
quite general applicability, and so it is plausible that a 
more sophisticated model might exhibit further bifurca- 
tions, thus producing different reduced periods and oscil- 
latory regimes. It may also be useful to bear in mind in 
this connection that three period halvings would result in 
11 years/2 3 ~ 1.3 years! We shall return to a more detailed 
study of the underlying dynamics as well as a quantita- 
tive study of different dynamo models elsewhere. We have 
chosen P r = 1 in order to obtain larger amplitude tor- 
sional oscillations near the surface. We have checked that 
fragmentation still occurs at smaller values of P r . 

Inevitably the uncertainties associated with the inver- 
sion of the helioseismic data so deep in the convection zone 
are quite large. Thus we believe that the mechanism dis- 
cussed here may, by demonstrating what modes of dynam- 
ical behaviour are theoretically possible, act as a concep- 
tual aid in interpreting current and further observations. 
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